// QtCore
#include <QtCore/QDir>
#include <QtCore/QFile>
#include <QtCore/QTextCodec>
#include <QtCore/QDataStream>
#include <QtCore/QStringBuilder>
#include <QtCore/QCoreApplication>
// Локальные
#include "settings.h"
#include "math/math.h"
#include "linal/los.h"
#include "linal/bcg.h"
#include "solver_tct.h"
#include "fem/trianglefem.h"
#include "geometry/2d/vector.h"

//================================================================================================================================
//================================================================================================================================
SolverTCt::SolverTCt() : m_NetSearcher(NULL)
{
    //****************************************************************************************************************************
    // Создание и очистка директории для вывода результатов.
    //****************************************************************************************************************************
    m_OutputPath = Settings::instance()->setting("OutputPath").toString();
    if(!QDir::isAbsolutePath(m_OutputPath))
        m_OutputPath = QDir::cleanPath(m_OutputPath).prepend(qApp->applicationDirPath() + "/");

    if(!QDir(m_OutputPath).exists())
        QDir(qApp->applicationDirPath()).mkpath(m_OutputPath);
    if(!QDir(m_OutputPath + "/txt").exists())
        QDir(qApp->applicationDirPath()).mkpath(m_OutputPath + "/txt");
    // TODO: Очистка директории.

    //****************************************************************************************************************************
    // Считывание параметров временной сетки.
    //****************************************************************************************************************************
    m_NumOfTimeLayers = Settings::instance()->setting("Time/TimeLayers").toInt();
    m_TimeStep = Settings::instance()->setting("Time/TimeStep").toDouble();

    //****************************************************************************************************************************
    // Считывание физических и геометрических параметров задачи. Перевод скорости из см/ч в м/с.
    //****************************************************************************************************************************
    m_TCrystal = Settings::instance()->setting("Physics/CrystallizationTemperature").toDouble();
    m_TMax = Settings::instance()->setting("Physics/MaxTemperature").toDouble();
    m_TRod = Settings::instance()->setting("Physics/RodTemperature").toDouble();

    m_TGrad = Settings::instance()->setting("Physics/TemperatureGradient").toDouble();
    m_TGrad *= 100.;    //К/см -> К/м

    m_V = Settings::instance()->setting("Physics/CrucibleSpeed").toDouble();
    m_V /= 100.;        //см/ч -> м/ч
    m_V /= 3600.;       //м/ч -> м/с

    //****************************************************************************************************************************
    // Считывание параметров решателей.
    //****************************************************************************************************************************
    m_MaxiterSLAE = Settings::instance()->setting("Math/SlaeIters").toInt();
    m_EpsSLAE = Settings::instance()->setting("Math/SlaeResidual").toDouble();

    //****************************************************************************************************************************
    // Инициализация решателей.
    //****************************************************************************************************************************
    m_Solvers.append(new linal::LOS(m_MaxiterSLAE, m_EpsSLAE));
    m_Solvers.append(new linal::BCG(m_MaxiterSLAE, m_EpsSLAE));

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    m_FrontHalfWidth = Settings::instance()->setting("Math/FrontHalfwidth").toDouble();
    m_MinCrystalHeight = Settings::instance()->setting("Math/MinCrystalHeight").toDouble();

    //****************************************************************************************************************************
    // Считывание сервисных настроек программы.
    //****************************************************************************************************************************
    m_OutputStep = Settings::instance()->setting("Solver/OtputStep").toInt();
}

//================================================================================================================================
//================================================================================================================================
bool SolverTCt::solveSLAE(const linal::CSRMatrix &A, const linal::Vector &f, linal::Vector &x)
{
    linal::Vector x_old = x;

    foreach(linal::SLAESolver* solver, m_Solvers)
    {
        x = x_old;
        if(solver->solveLU(A,f,x) < m_EpsSLAE)
            return true;
    }

    foreach(linal::SLAESolver* solver, m_Solvers)
    {
        x = x_old;
        if(solver->solveDi(A,f,x) < m_EpsSLAE)
            return true;
    }

    foreach(linal::SLAESolver* solver, m_Solvers)
    {
        x = x_old;
        if(solver->solve(A,f,x) < m_EpsSLAE)
            return true;
    }

    qWarning("SLAE wasn't solved!");
    return false;
}

//================================================================================================================================
//================================================================================================================================
void SolverTCt::solve()
{
    QTextStream cout(stdout);
    cout.setCodec(QTextCodec::codecForName("IBM866"));

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    int resume_from = Settings::instance()->setting("Solver/ResumeFrom").toInt();

    if(resume_from != -1)
        readResumeInfo(resume_from);

    int first_layer = (resume_from == -1 ? 1 : resume_from + 1);
    for(int cur_layer = first_layer; cur_layer < m_NumOfTimeLayers; cur_layer++)
    {
        //Решаем уравнение для T
        m_A.zeroize();
        m_F.zeroize();
        formTSLAE(cur_layer);
        solveSLAE(m_A, m_F, m_T);

        recreateNet();

        m_A.zeroize();
        m_F.zeroize();
        formCSLAE(cur_layer);
        solveSLAE(m_A, m_F, m_C);

        m_TOldTime = m_T;
        m_COldTime = m_C;

        m_OldFront.clear();
        for(int i = 0; i < m_Front.size(); ++i)
            m_OldFront.append(m_Net.point(m_Front[i]));

        if(cur_layer%m_OutputStep == 0)
            writeResult(cur_layer);
        if(cur_layer%50 == 0)
            cout << cur_layer << endl;
    }
}

//================================================================================================================================
//================================================================================================================================
void SolverTCt::tBC(int time_layer)
{
    for(int i = 0; i < m_TopTBC.size(); i++)
    {
        m_A.makeIdentity(m_TopTBC[i]);

        double time = time_layer*m_TimeStep;
        switch(m_TypeTBC[i])
        {
        case 0:
            m_F[m_TopTBC[i]] = m_TMax - m_V*m_TGrad*time;
            break;
        case 1:
            if(m_Net.point(m_TopTBC[i]).y() >= time*m_V)
                m_F[m_TopTBC[i]] = m_TMax;
            else
                m_F[m_TopTBC[i]] = m_TGrad*m_Net.point(m_TopTBC[i]).y() + m_TMax - m_V*m_TGrad*time;
            break;
        case 2:
            m_F[m_TopTBC[i]] = m_TRod;
            break;
        default:
            qFatal("Unknown BC type");
        }
    }
}

//================================================================================================================================
//================================================================================================================================
void SolverTCt::cBC(int time_layer)
{
    Q_UNUSED(time_layer);

    if(m_State == SolverTCt::JustStarted || m_State == SolverTCt::InProgress)
    {
        QVector<double> speed = frontSpeed();

        for(int i = 0; i < m_Front.size() - 1; i++)
        {
            double dr = m_Net.point(m_Front[i + 1]).x() - m_Net.point(m_Front[i]).x();
            double dz = m_Net.point(m_Front[i + 1]).y() - m_Net.point(m_Front[i]).y();
            double h = sqrt(dr*dr + dz*dz);
            double r1 = m_Net.point(m_Front[i]).x();
            double r2 = m_Net.point(m_Front[i + 1]).x();
            double v1 = speed[i];
            double v2 = speed[i + 1];
            double v = (v1 + v2)/2.;

            m_A(m_Front[i], m_Front[i]) += (r1/4. + r2/12.)*h*v*(1 - m_Materials[5].property(Material::SegregationFactor));
            m_A(m_Front[i + 1], m_Front[i + 1]) += (r1/12. + r2/4.)*h*v*(1 - m_Materials[5].property(Material::SegregationFactor));

            m_A(m_Front[i], m_Front[i + 1]) += (r1/12. + r2/12.)*h*v*(1 - m_Materials[5].property(Material::SegregationFactor));
            m_A(m_Front[i + 1], m_Front[i]) += (r1/12. + r2/12.)*h*v*(1 - m_Materials[5].property(Material::SegregationFactor));
        }
    }
}

//================================================================================================================================
//================================================================================================================================
void SolverTCt::formTSLAE(int time_layer)
{
    double density;
    double heat_capacity;
    double heat_conductivity;
    double D[3][3];             // Матрица, используемая для определения базисных функций
    double C[3][3][3];          // Локальная матрица массы
    double detD;                // Определитель матрицы D
    double abs_detD;            // Модуль определетиля матрицы D
    double ht;                  // Шаг по времени
    double r[3];                // r-координаты вершин КЭ
    double z[3];                // z-координаты вершин КЭ
    int top[3];                 // Индексы вершин КЭ

    ht = m_TimeStep;            // Задел на будущее, когда сетка может быть неравномерной.

    //****************************************************************************************************************************
    //Собираем конечноэлементную СЛАУ
    //****************************************************************************************************************************
    for(int cur_triangle = 0; cur_triangle < m_Net.numOfElements(); cur_triangle++)
    {
        density = m_Materials[m_Net.material(cur_triangle)].property(Material::Density);
        heat_capacity = m_Materials[m_Net.material(cur_triangle)].property(Material::HeatCapacity);
        heat_conductivity = m_Materials[m_Net.material(cur_triangle)].property(Material::HeatConductivity);

        for(int i = 0; i < 3; i++)
        {
            r[i] = m_Net.point(cur_triangle, i).x();
            z[i] = m_Net.point(cur_triangle, i).y();
            top[i] = m_Net.pointIndex(cur_triangle, i);
        }

        detD = (r[1] - r[0])*(z[2] - z[0]) - (r[2] - r[0])*(z[1] - z[0]);
        abs_detD = qAbs(detD);

        D[0][0] = (r[1]*z[2] - r[2]*z[1])/detD;
        D[0][1] = (z[1] - z[2])/detD;
        D[0][2] = (r[2] - r[1])/detD;
        D[1][0] = (r[2]*z[0] - r[0]*z[2])/detD;
        D[1][1] = (z[2] - z[0])/detD;
        D[1][2] = (r[0] - r[2])/detD;
        D[2][0] = (r[0]*z[1] - r[1]*z[0])/detD;
        D[2][1] = (z[0] - z[1])/detD;
        D[2][2] = (r[1] - r[0])/detD;

        C[0][0][0] = abs_detD/20.;      C[0][0][1] = abs_detD/60.;      C[0][0][2] = abs_detD/60.;
        C[0][1][0] = abs_detD/60.;      C[0][1][1] = abs_detD/60.;      C[0][1][2] = abs_detD/120.;
        C[0][2][0] = abs_detD/60.;      C[0][2][1] = abs_detD/120.;     C[0][2][2] = abs_detD/60.;

        C[1][0][0] = abs_detD/60.;      C[1][0][1] = abs_detD/60.;      C[1][0][2] = abs_detD/120.;
        C[1][1][0] = abs_detD/60.;      C[1][1][1] = abs_detD/20.;      C[1][1][2] = abs_detD/60.;
        C[1][2][0] = abs_detD/120.;     C[1][2][1] = abs_detD/60.;      C[1][2][2] = abs_detD/60.;

        C[2][0][0] = abs_detD/60.;      C[2][0][1] = abs_detD/120.;     C[2][0][2] = abs_detD/60.;
        C[2][1][0] = abs_detD/120.;     C[2][1][1] = abs_detD/60.;      C[2][1][2] = abs_detD/60.;
        C[2][2][0] = abs_detD/60.;      C[2][2][1] = abs_detD/60.;      C[2][2][2] = abs_detD/20.;

        for(int i = 0; i < 3; i++)
        {
            for(int j = 0; j < 3; j++)
            {
                double addition(0);
                addition += (r[0] + r[1] + r[2])*(D[i][1]*D[j][1] + D[i][2]*D[j][2])*heat_conductivity*abs_detD/6.;
                addition += (C[i][j][0]*r[0] + C[i][j][1]*r[1] + C[i][j][2]*r[2])*heat_capacity*density/ht;
                m_A(top[i], top[j]) += addition;
            }
        }

        for(int i = 0; i < 3; i++)
            for(int j = 0; j < 3; j++)
                for(int k = 0; k < 3; k++)
                    m_F[top[i]] += m_TOldTime[top[j]]*r[k]*C[i][j][k]*heat_capacity*density/ht;
    }

    //****************************************************************************************************************************
    //Учитываем краевые условия.
    //****************************************************************************************************************************
    tBC(time_layer);
}

//================================================================================================================================
//================================================================================================================================
void SolverTCt::formCSLAE(int time_layer)
{
    double diffusivity;
    double D[3][3];         // Матрица, используемая для определения базисных функций
    double C[3][3][3];      // Локальная матрица массы
    double detD;            // Определитель матрицы D
    double abs_detD;        // Модуль определетиля матрицы D
    double ht;              // Шаг по времени
    double r[3];            // r-координаты вершин КЭ
    double z[3];            // z-координаты вершин КЭ
    int top[3];             // Индексы вершин КЭ

    ht = m_TimeStep;        // Задел на будущее, когда сетка может быть неравномерной.

    //****************************************************************************************************************************
    //Собираем конечноэлементную СЛАУ
    //****************************************************************************************************************************
    for(int cur_triangle = 0; cur_triangle < m_Net.numOfElements(); cur_triangle++)
    {
        if(m_Materials[m_Net.material(cur_triangle)].state() == Material::Liquid)
        {
            diffusivity = m_Materials[5].property(Material::Diffusivity);

            for(int i = 0; i < 3; i++)
            {
                r[i] = m_Net.point(cur_triangle, i).x();
                z[i] = m_Net.point(cur_triangle, i).y();
                top[i] = m_Net.pointIndex(cur_triangle, i);
            }

            detD = (r[1] - r[0])*(z[2] - z[0]) - (r[2] - r[0])*(z[1] - z[0]);
            abs_detD = qAbs(detD);

            D[0][0] = (r[1]*z[2] - r[2]*z[1])/detD;
            D[0][1] = (z[1] - z[2])/detD;
            D[0][2] = (r[2] - r[1])/detD;
            D[1][0] = (r[2]*z[0] - r[0]*z[2])/detD;
            D[1][1] = (z[2] - z[0])/detD;
            D[1][2] = (r[0] - r[2])/detD;
            D[2][0] = (r[0]*z[1] - r[1]*z[0])/detD;
            D[2][1] = (z[0] - z[1])/detD;
            D[2][2] = (r[1] - r[0])/detD;

            C[0][0][0] = abs_detD/20.;      C[0][0][1] = abs_detD/60.;      C[0][0][2] = abs_detD/60.;
            C[0][1][0] = abs_detD/60.;      C[0][1][1] = abs_detD/60.;      C[0][1][2] = abs_detD/120.;
            C[0][2][0] = abs_detD/60.;      C[0][2][1] = abs_detD/120.;     C[0][2][2] = abs_detD/60.;

            C[1][0][0] = abs_detD/60.;      C[1][0][1] = abs_detD/60.;      C[1][0][2] = abs_detD/120.;
            C[1][1][0] = abs_detD/60.;      C[1][1][1] = abs_detD/20.;      C[1][1][2] = abs_detD/60.;
            C[1][2][0] = abs_detD/120.;     C[1][2][1] = abs_detD/60.;      C[1][2][2] = abs_detD/60.;

            C[2][0][0] = abs_detD/60.;      C[2][0][1] = abs_detD/120.;     C[2][0][2] = abs_detD/60.;
            C[2][1][0] = abs_detD/120.;     C[2][1][1] = abs_detD/60.;      C[2][1][2] = abs_detD/60.;
            C[2][2][0] = abs_detD/60.;      C[2][2][1] = abs_detD/60.;      C[2][2][2] = abs_detD/20.;

            for(int i = 0; i < 3; i++)
            {
                for(int j = 0; j < 3; j++)
                {
                    double addition(0);
                    addition += (r[0] + r[1] + r[2])*(D[i][1]*D[j][1] + D[i][2]*D[j][2])*diffusivity*abs_detD/6.;
                    addition += (C[i][j][0]*r[0] + C[i][j][1]*r[1] + C[i][j][2]*r[2])/ht;
                    m_A(top[i], top[j]) += addition;
                }
            }

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    for(int k = 0; k < 3; k++)
                        m_F[top[i]] += m_COldTime[top[j]]*r[k]*C[i][j][k]/ht;
        }
    }

    //*******************************************************************************************************
    //Учитываем краевые условия
    //*******************************************************************************************************
    cBC(time_layer);

    //*******************************************************************************************************
    //Заполняем матрицу для нерасплава
    //*******************************************************************************************************
    for(int i = 0; i < m_A.dimension(); ++i)
        if(qAbs(m_A(i,i)) < DBL_EPSILON)
        {
            m_A(i, i) = 1.;
            m_F[i] = 0.;
        }
}

//================================================================================================================================
//Вывод результата в файл.
//================================================================================================================================
void SolverTCt::writeResult(int time_layer)
{
    QFile file_bin;
    QFile file_txt;
    QDataStream ostream_bin(&file_bin);
    QTextStream ostream_txt(&file_txt);
    ostream_bin.setVersion(QDataStream::Qt_4_7);
    ostream_txt.setRealNumberPrecision(12);

    // Количество десятичных разрядов в числе, равном количеству временных слоев.
    // Не все программы способны "правильно" упорядочивать файлы, в именах которых встречаются числа с разным количеством разрядов.
    // Так, например, одни программы считают, что файл 10.dat должен распологаться раньше файла 100.dat, в то время как другие
    // программы наоборот располагают файл 100.dat перед файлом 10.dat. Единственный способ заставить все программы одинаково
    // работать с подобными файлами - выравнять количество разрядов в числах, присутствующих в именах файлов, добавляя при
    // необходимости незначащие нули перед числом. Т.е. файл 010.dat любой программой будет расположен перед файлом 100.dat.
    // В данном приложении, в имени файлов указываются номер временного слоя и, сл-но, для того, чтобы знать нужно ли добавлять
    // незначащие нули, и если да, то сколько, необходимо знать число разрядов в максимально возможном номере временного слоя.
    int num_of_orders = QString::number(m_NumOfTimeLayers).count();

    //****************************************************************************************************************************
    //Выводим конечноэлементую сетку.
    //****************************************************************************************************************************
    QString filename_nvtr = QString("/%1nvtr.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);
    QString filename_rz = QString("/%1rz.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);
    m_Net.writeToFiles(filename_rz, filename_nvtr);

    file_txt.setFileName(QString("/txt/%1nvtr.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_txt.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfElements(); ++i)
    {
        for(int j = 0; j < 3; ++j)
            ostream_txt << m_Net.pointIndex(i, j) << " ";
        ostream_txt << m_Net.material(i) << endl;
    }
    file_txt.close();

    file_txt.setFileName(QString("/txt/%1rz.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_txt.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfPoints(); ++i)
        ostream_txt << m_Net.point(i) << endl;
    file_txt.close();

    //****************************************************************************************************************************
    //Выводим значения температуры.
    //****************************************************************************************************************************
    file_bin.setFileName(QString("/%1T.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_bin.open(QFile::WriteOnly);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream_bin << m_Net.point(i).x() << m_Net.point(i).y() << m_T[i];
    file_bin.close();

    file_txt.setFileName(QString("/txt/%1T.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_txt.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream_txt << m_Net.point(i) << " " << m_T[i] << endl;
    file_txt.close();

    //****************************************************************************************************************************
    //Выводим значения концентрации.
    //****************************************************************************************************************************
    file_bin.setFileName(QString("/%1C.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_bin.open(QFile::WriteOnly);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream_bin << m_Net.point(i).x() << m_Net.point(i).y() << m_C[i];
    file_bin.close();

    file_txt.setFileName(QString("/txt/%1C.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_txt.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream_txt << m_Net.point(i) << " " << m_C[i] << endl;
    file_txt.close();

    //************************************************************************************************************************
    //выводим информацию для восстановления решения.
    //************************************************************************************************************************
    file_bin.setFileName(QString("/%1bc.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_bin.open(QIODevice::WriteOnly);
    ostream_bin << m_TopTBC.size();
    for(int i = 0; i < m_TopTBC.size(); i++)
        ostream_bin << m_TopTBC[i] << m_TypeTBC[i];
    file_bin.close();

    file_bin.setFileName(QString("/%1front.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_bin.open(QIODevice::WriteOnly);
    ostream_bin << m_Front.size();
    for(int i = 0; i < m_Front.size(); i++)
        ostream_bin << m_Net.point(m_Front[i]).x() << m_Net.point(m_Front[i]).y();
    file_bin.close();

    file_bin.setFileName(QString("/%1state.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_bin.open(QIODevice::WriteOnly);
    ostream_bin << m_State;
    file_bin.close();
}

//================================================================================================================================
//Инициализация структур данных, необходимая для продолжения расчетов с указанного места.
//================================================================================================================================
void SolverTCt::readResumeInfo(int time_layer)
{
    // Количество десятичных разрядов в числе, равном количеству временных слоев.
    // Не все программы способны "правильно" упорядочивать файлы, в именах которых встречаются числа с разным количеством разрядов.
    // Так, например, одни программы считают, что файл 10.dat должен распологаться раньше файла 100.dat, в то время как другие
    // программы наоборот располагают файл 100.dat перед файлом 10.dat. Единственный способ заставить все программы одинаково
    // работать с подобными файлами - выравнять количество разрядов в числах, присутствующих в именах файлов, добавляя при
    // необходимости незначащие нули перед числом. Т.е. файл 010.dat любой программой будет расположен перед файлом 100.dat.
    // В данном приложении, в имени файлов указываются номер временного слоя и, сл-но, для того, чтобы знать нужно ли добавлять
    // незначащие нули, и если да, то сколько, необходимо знать число разрядов в максимально возможном номере временного слоя.
    int num_of_orders = QString::number(m_NumOfTimeLayers).count();

    QString filename_rz = QString("/%1rz.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);
    QString filename_nvtr = QString("/%1nvtr.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);

    QString filename_t = QString("/%1T.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);
    QString filename_c = QString("/%1C.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);

    QString filename_bc = QString("/%1bc.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);
    QString filename_front = QString("/%1front.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);
    QString filename_state = QString("/%1state.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    m_Net.readFromFiles(filename_rz, filename_nvtr);
    m_NetSearcher->analyzeNet(m_Net);

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    int data_size;
    int bc_top, bc_type;

    QFile file(filename_bc);
    QDataStream istream(&file);
    istream.setVersion(QDataStream::Qt_4_7);
    file.open(QIODevice::ReadOnly);

    istream >> data_size;
    m_TopTBC.clear();
    m_TypeTBC.clear();
    for(int i = 0; i < data_size; i++)
    {
        istream >> bc_top >> bc_type;
        m_TopTBC.append(bc_top);
        m_TypeTBC.append(bc_type);
    }

    file.close();

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    double r,z;

    file.setFileName(filename_front);
    file.open(QIODevice::ReadOnly);

    istream >> data_size;
    m_OldFront.clear();
    for(int i = 0; i < data_size; i++)
    {
        istream >> r >> z;
        m_OldFront.append(geometry::_2d::Point(r, z));
    }

    file.close();

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    int state;

    file.setFileName(filename_state);
    file.open(QIODevice::ReadOnly);
    istream >> state;
    switch(state)
    {
    case 0: m_State = SolverTCt::NotStarted; break;
    case 1: m_State = SolverTCt::JustStarted; break;
    case 2: m_State = SolverTCt::InProgress; break;
    case 3: m_State = SolverTCt::Finished; break;
    default: qFatal("Unknown crystallization state");
    }
    file.close();

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    m_T.readFromFile(filename_t);
    m_TOldTime = m_T;

    m_C.readFromFile(filename_c);
    m_COldTime = m_C;

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    m_F.resize(m_Net.numOfPoints());

    QVector<int> igl, jgl, igu, ngu;
    fem::triangles::formMatrixProfile(m_Net, igl, jgl, igu, ngu);
    m_A = linal::CSRMatrix(m_Net.numOfPoints(), igl, jgl, igu, ngu);
}
